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The  inverse  design  problem  technique  presented  in  this  paper  is  intended  for  optimizing  the  shape  of  the 
gas  channel  at  the  cathode  side  in  a  proton  exchange  membrane  fuel  cell  (PEMFC).  The  technique  uses  the 
desired  current  densities  located  on  a  carbon  plate  near  the  outlet  of  the  channel  at  the  cathode  side  as  a 
starting  point.  The  desired  current  density  distributions  can  be  obtained  by  modifying  the  current  density 
distributions  of  the  existing  PEMFC  with  rectangular  gas  channels.  The  geometry  of  the  redesigned  gas 
channel  is  generated  using  a  B-spline  curve  method,  which  enables  the  shape  of  the  fuel  channel  to  be 
completely  specified  using  only  a  small  number  of  control  points,  thus  applying  the  technique  of  parameter 
estimation  for  the  inverse  design  problem.  Results  show  that  by  utilizing  the  redesigned  optimal  gas 
channel,  the  total  current  of  the  PEMFC  can  be  increased,  and  at  the  same  time  the  phenomenon  of 
saturated  water  accumulation  in  the  channel  can  be  greatly  reduced. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

In  order  to  adapt  to  environmental  change  as  well  as  the 
world’s  current  energy  crisis,  in  the  past  several  years  researchers 
around  the  world  have  been  focusing  on  developing  renewable 
resources  including  geothermal  energy,  solar  energy,  wind  energy, 
tidal  energy  and  hydrogen  energy  to  substitute  for  conventional 
resources  such  as  coal,  oil,  and  natural  gas.  Fuel  cells  have  sev¬ 
eral  advantages  over  conventional  power  generation  appliances, 
and  can  thus  be  considered  major  power-generating  devices  in 
the  future.  In  a  fuel  cell,  chemical  energy  is  converted  to  electrical 
energy  directly,  with  little  heat  production.  This  implies  that  the 
efficiency  of  a  fuel  cell  is  not  limited  by  the  Carnot  cycle  efficiency. 

In  general,  fuel  cells  produce  less  toxic  by-products  than  con¬ 
ventional  appliances,  thus  making  them  environment-friendly. 
Moreover,  fuel  cells  are  very  quiet  during  operation  and  can  be  of 
various  sizes;  these  features  will  make  them  ideal  for  installation 
in  any  place.  The  application  range  of  fuel  cells  is  very  wide,  from 
a  soldier’s  equipment  on  the  battlefield  to  huge  power  generators 
for  residential  and  industrial  applications.  Most  fuel  cells  operate 
at  low  temperatures  (around  70  °C),  with  the  exception  of  solid 
oxide  fuel  cells,  and  this,  added  to  the  fact  that  they  do  not  have 
any  moving  parts,  will  make  them  easy  to  design  from  a  materials 
standpoint. 
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Significant  progress  in  the  development  of  fuel  cell  technol¬ 
ogy  has  been  achieved  by  an  increasing  number  of  experimental 
and  theoretical  studies  [1-6].  One  of  the  most  promising  types 
of  fuel  cells,  the  proton  exchange  membrane  fuel  cell  (PEMFC), 
is  currently  being  aggressively  researched  due  to  its  significant 
advantages  in  portable  electronic  applications  over  conventional 
battery  systems.  It  also  meets  more  stringent  emissions  stan¬ 
dards,  because  far  fewer  pollutants  are  produced  by  it  than 
by  internal  combustion  engines.  Common  flow-field  designs  for 
this  fuel  cell  include  the  conventional  flow  fields,  i.e.,  parallel 
flow  field,  Z-type  flow  field,  and  serpentine  flow  field.  Due  to 
the  fact  that  flow-field  design  in  the  fuel  cell  is  important  in 
determining  the  fuel  cell’s  performance,  many  researchers  have 
examined  this  issue  in  recent  years.  However,  for  simplicity,  in 
the  present  study  only  the  case  of  the  parallel  flow  field  is  exam¬ 
ined. 

In  the  PEMFC,  the  maximum  current  densities  will  appear  at 
the  flow  channel  inlets  and  decrease  along  the  channels  to  the  out¬ 
lets.  This  is  due  to  the  fact  that  the  oxygen  is  gradually  consumed 
and  liquid  water  is  accumulated  due  to  the  electrochemical  reaction 
along  the  flow  channels.  For  this  reason,  if  the  large  amount  of  liquid 
water  accumulation  can  be  removed  efficiently,  the  current  densi¬ 
ties  along  the  channels  will  be  increased  and  the  performance  of  the 
PEMFC  must  therefore  be  improved  significantly.  Many  efforts  have 
been  devoted  to  optimizing  the  flow  field  design  in  order  to  improve 
the  cell  performance  [7-12].  This  can  also  be  done  by  redesigning 
the  shape  of  the  gas  channel,  applying  the  technique  of  the  inverse 
design  problem  [13-15]. 
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The  direct  PEMFC  problems  are  concerned  with  the  determina¬ 
tion  of  the  current  densities,  saturated  water  and  pressure  drops  in 
the  cell  when  the  operational  conditions,  system  parameters  and 
the  shape  of  gas  channels  are  all  specified.  In  contrast,  the  inverse 
geometry  design  problem  for  the  PEMFC  considered  in  this  paper 
involves  the  determination  of  the  optimal  shape  of  gas  channels 
based  on  the  requirement  of  desired  current  density  on  the  upper 
carbon  plate  of  the  redesigned  gas  channel  at  the  cathode. 

Recently,  Cheng  et  al.  [16]  applied  the  simplified  conjugate  gra¬ 
dient  method  (SCGM)  in  the  PEMFC  for  the  optimization  of  the 
geometric  parameters  of  the  fuel  cells  and  obtained  good  esti¬ 
mations.  Yan  et  al.  [17]  proposed  a  parallel  flow  channel  design 
with  tapered  heights  or  widths  to  improve  reactant  utilization  effi¬ 
ciency  for  PEMFCs.  In  their  study,  the  tapered  channel  designs  are 
used  and  the  flow  area  contraction  along  the  flow  channel  leads  to 
increased  reactant  velocities,  thereby  enhancing  the  reactant  trans¬ 
port  through  the  porous  layers,  reactant  utilization  and  liquid  water 
removal.  However,  one  should  note  that  these  tapered  channels 
are  designed  manually  and  are  not  based  on  the  algorithm  of  the 
inverse  design  problem. 

With  this  in  mind,  the  technique  of  the  inverse  design  problem 
should  be  used  to  design  the  optimal  gas  channel  geometry  in  accor¬ 
dance  with  the  desired  current  density  distributions  of  the  cell.  The 
desired  current  density  distributions  can  be  obtained  by  modifying 
the  existing  current  density  distributions  in  the  PEMFC.  A  B-spline 
curve  is  utilized  to  simulate  the  shape  of  the  gas  channel,  and  there¬ 
fore  the  shape  parameters  in  this  formulation  are  the  coordinates 
of  each  control  point,  i.e.,  this  limited  number  of  control  points 
becomes  the  set  of  parameters  controlling  the  gas  channel  geom¬ 
etry.  In  the  present  study,  the  task  is  to  redesign  the  gas  channel 
geometry  for  parallel  flow  arrangement,  based  on  both  the  orig¬ 
inal  gas  channel  geometry  and  the  new  desired  current  density 
distributions.  The  originality  of  this  study  lies  in  that  the  present 
algorithm  can  design  the  optimal  shape  of  gas  channel  that  satisfies 
the  desired  current  density. 

For  the  present  geometry  design  problem,  due  to  its  inher¬ 
ent  nature,  a  complete  regeneration  of  the  mesh  is  required  as 
the  geometry  evolves.  Moreover,  the  continuous  evolution  of  the 
geometry  itself  poses  certain  difficulties  in  arriving  at  analytical  or 
numerical  solutions.  For  this  reason,  it  is  necessary  to  use  an  effi¬ 
cient  technique  that  can  handle  a  problem  with  irregular  surface 
geometry,  especially  in  this  3D  application. 

This  paper  addresses  the  development  of  an  efficient  method 
for  parameter  estimation,  i.e.,  the  Levenberg-Marquardt  algorithm, 
in  estimating  the  optimal  gas  channel  geometry  that  satisfies  the 
desired  current  density  distributions.  The  Levenberg-Marquardt 
method  has  proven  to  be  a  powerful  algorithm  in  inverse  design 
calculations.  This  inverse  design  method  had  been  applied  to  pre¬ 
dict  the  form  of  a  ship’s  hull  in  accordance  with  the  desired  hull 
pressure  distribution  by  Huang  et  al.  [13].  Subsequently,  Chen  and 
Huang  [14]  applied  it  to  predict  an  unknown  hull  form  based  on  the 
preferable  wake  distribution  in  the  propeller  disk  plane.  Chen  et  al. 
[15]  further  applied  it  to  the  aspect  of  optimal  design  for  a  bulbous 
bow. 

In  Section  2  of  this  paper,  the  algorithm  used  to  calculate  the  cell 
current  density  distribution,  i.e.,  the  direct  problem,  is  explained. 
The  method  of  B-spline  curves  is  described  in  Section  3.  The  inverse 
design  problem,  involving  the  definition  of  a  cost  function  and  the 
Levenberg-Marquardt  algorithm,  is  addressed  in  Section  4.  Finally, 
a  computational  procedure  is  summarized  in  Section  5. 

2.  The  direct  problem 

Fuel  cells  constitute  an  important  alternative  to  power  gener¬ 
ation  by  engines  or  turbines.  With  an  increasing  trend  towards 


reduced  emissions  and  higher  efficiency,  fuel  cells  are  emerging 
rapidly  as  the  technology  of  choice  for  power  generation  in  the  new 
century. 

The  PEMFC  consists  of  an  anode  flow  channel,  an  anode  dif¬ 
fuser  layer,  an  anode  catalyst  layer,  a  proton  exchange  membrane,  a 
cathode  catalyst  layer,  a  cathode  diffuser  layer,  and  a  cathode  flow 
channel.  To  simplify  the  analysis,  it  is  assumed  that  the  gas  mixtures 
are  ideal  gases;  the  flow  field  is  steady,  laminar  and  incompressible; 
the  cell  is  isothermal;  the  electrochemical  reaction  only  takes  place 
in  the  catalyst  layer;  and  the  cell  structures  in  the  diffuser  layers, 
the  catalyst  layers,  and  the  proton  exchange  membrane  are  made 
of  isotropic  and  homogeneous  porous  material. 

The  gas  channels  and  the  gas  diffusion  layers  inside  the  fuel 
cells  are  designed  to  allow  the  reactant  gases  to  diffuse  uniformly 
into  the  catalyst  layers,  and  at  the  same  time,  to  prevent  excessive 
increase  in  electrical  resistance  against  the  released  electrons  from 
the  catalyst  layers.  The  porous  layers  include  the  gas  diffusion  layers 
and  the  catalyst  layers,  in  which  the  mass,  momentum,  and  species 
equations  for  the  flows  should  be  derived  based  on  the  non-Darcy 
law. 

The  mathematical  models  for  the  present  3D  PEMFC  involving 
four  groups  of  differential  equations,  namely,  (a)  mass,  momentum, 
and  species  transport  phenomena  in  the  gas  channels  (GC)  and  the 
porous  layers— i.e.,  gas  diffusion  layer  (GDL)  and  catalyst  layer  (CL); 
(b)  electrochemical  reactions  in  the  catalyst  layers;  (c)  ionic  con¬ 
duction  taking  place  in  the  catalyst  layers  and  the  proton  exchange 
membrane;  (d)  electronic  conduction  in  the  porous  layers  and  the 
carbon  plates  (CP).  These  governing  equations  are  described  below 
[16,17]. 

2.1.  Mass,  momentum,  and  species  transport  equations 

The  conservation  equations  for  mass,  momentum,  and  species 
can  be  used  to  model  the  gas  flows  in  the  gas  channels,  gas  diffusion 
layer  and  catalyst  layer  for  the  PEMFC.  These  equations  are  given 
below  for  mass,  momentum  and  species  equations,  respectively: 

V-(£pL0  =  O,  (1) 

-  p2ulJ 

V  •  (spUU)  =  -sVP  +  V  •  ( sx )  +  — C— ,  (2) 

K 

V  •  ( epUZn )  =  V  •  (DfVZn )  +  S„ .  (3 ) 

Here  U  denotes  the  gas  velocity  vector,  e  and  k  are  the  porosity 
and  permeability  of  the  porous  layers,  respectively,  and  r  is  the  fluid 
stress  tensor.  The  fuel  gases  at  the  anode  side  and  cathode  side  are 
H2  and  02,  respectively.  Zn  represents  the  species  mass  fraction  of 
species  n.  It  should  be  noted  that  e  =  1  and  /<  =  oo  imply  pure  fluid 
flows  in  the  gas  channels. 

The  effective  mass  diffusivity  D£ff  of  species  n  can  be  obtained 
by  using  Bruggemann’s  equation,  and  is  calculated  as 

Df  =  eeDn ,  (4) 

where  0  is  the  Bruggemann  coefficient  of  the  porous  layer  and  Dn  is 
the  mass  diffusivity  of  species  n.  In  Eq.  (3),  Sn  represents  the  source 
term  of  species  n,  which  is  due  to  the  electrochemical  reactions  in 
the  catalyst  layers. 

2.2.  Electrochemical  reaction  equations 

During  the  reaction,  H2  and  02  are  consumed  in  the  catalyst  lay¬ 
ers  of  the  anode  and  cathode,  respectively,  and  H20  is  the  product 
of  the  reaction  at  the  cathode.  The  source  terms  for  different  species 
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are  defined  as  follows: 


of  water  activity  in  the  membrane  can  be  obtained  as 


c  =ZU 
h2  2 F  ’ 


Sq2  = 


k_ 

4F’ 


(5) 

(6) 


(7) 


where  F  is  the  Faraday  constant,  while  /a  and  Jc  are  the  transfer 
current  densities  at  the  anode  and  the  cathode,  respectively.  These 
values  can  be  obtained  based  on  the  Butler-Volmer  condition: 


where  /aref  and  /cref  are  the  reference  transfer  current  densities, 
and  ra  and  rc  the  concentration  parameters  at  the  anode  and  the 
cathode,  respectively.  The  terms  aa,a  and  aa<c  are  the  anodic  and 
the  cathodic  charge  transfer  coefficients  for  the  anode  reaction, 
while  ac,a  and  ac,c  are  the  anodic  and  the  cathodic  charge  trans¬ 
fer  coefficients  for  the  cathode  reaction.  Zh2  ref  and  Zq2  ref  are  the 
reference  concentrations  at  which  the  exchange  current  densities 
are  obtained. 


2.3.  Ionic  conduction  equations 

The  ionic  conduction  occurs  in  the  catalyst  layers  and  the  proton 
exchange  membrane,  since  protons  of  hydrogen  (H+)  are  gener¬ 
ated  in  the  anode  catalyst  layer  and  then  move  through  the  proton 
exchange  membrane  toward  the  cathode  catalyst  layer  to  react  with 
oxygen  and  the  returning  electrons.  The  ionic  conduction  equa¬ 
tion  for  determining  the  phase  potential,  cp,  in  the  catalyst  layers 
is  obtained  as 

-v.(rcfv^a)  =  fa,  (10) 

where  tpc i  and  r^[f  are  the  phase  potential  and  effective  ionic 
conductivity  of  the  catalyst  layers,  respectively.  The  effective  ionic 
conductivity,  7^f,  is  evaluated  by  using  the  following  equation: 

rcf  =  *&Lxn,a.  (11) 

The  effective  ionic  conductivity  of  the  catalyst  layer,  7^f, 
increases  with  the  porosity  of  solid  catalyst  particles  in  the  cata¬ 
lyst  layer,  £s  Cl.  The  ionic  conduction  equation  for  the  membrane 
can  be  expressed  as 


v-(rmv<?m)  =  o, 


(12) 


where  fm  is  the  ionic  conductivity  and  (pm  is  the  phase  potential 
of  the  membrane.  Springer  et  al.  [18]  suggested  that  the  expression 
for  ionic  conductivity  of  the  membrane  can  be  expressed  as 


—  T^m  ref  x  ^Xp 


,268*(353 


(13) 


where  the  reference  ionic  conductivity  rm  re f  is  given  as 


rm?ref  =  0.005239A  -  0.00326.  (14) 

Based  on  the  experimental  data,  the  water  content  of  the  mem¬ 
brane,  A,  depends  on  the  dimensionless  water  activity  in  the 
membrane,  la\ which  is  defined  as  the  gas  partial  pressure  of  water 
at  equilibrium  with  the  membrane  divided  by  the  saturation  pres¬ 
sure  of  water  vapor  at  the  operating  temperature.  The  expression 


a  = 


Xu2oP 

p^T' 


(15) 


The  saturation  pressure  of  water  vapor,  Psat,  can  be  computed 


by 

Psat  =  10 


=  Q-2.1794+0.02953(r— 273.15)— 9. 1837(T— 273. 15)2  xlO-5 +1.4454(7—273. 15)3  xlO-7 


with  T  in  K. 

The  relationship  between  A  and  a  can  be  obtained  as 

f  0.043  +  17.81a  -  39.85a2  +  36.0a3;  0  <  a  <  1 

A  =  <  14.0  +  1.4(a  - 1);  l<a<3  .  (17) 

[  16.8;  3  <  a 


Finally,  based  on  the  solution  for  the  phase  potential  tp,  the  dis¬ 
tribution  of  local  ionic  current  density,  /,  may  be  further  determined 
by  the  following  expression: 

l  =  -TV(p.  (18) 


2.4.  Electronic  conduction  equations 


The  carbon  plate,  the  GDL,  and  the  catalyst  layers  all  serve  as  con¬ 
ductors  for  electric  current.  In  the  carbon  plate  and  the  GDL,  there 
is  no  electrochemical  reaction  taking  place.  The  electrochemical 
reactions  only  take  place  in  the  catalyst  layers.  Due  to  the  electro¬ 
chemical  reactions  activated  in  the  catalyst  layers,  there  exists  a 
source  term  in  the  electronic  conduction  equation.  Therefore,  the 
electronic  conduction  equation  for  these  three  components  can  be 
expressed  in  the  following  general  form: 

V.(o fV0)  =  fd.  (19) 

Here,  ‘d’  represents  three  different  layers,  i.e.,  the  carbon  plate, 
the  GDL,  and  the  catalyst  layer,  respectively.  The  term  0  denotes  the 
electronic  potential  and  <r^ff  is  the  effective  electronic  conductivity. 
The  effective  electronic  conductivity  for  the  carbon  plate  (CP)  can 
be  taken  as  a  constant,  while  the  effective  electronic  conductivities 
for  the  porous  gas  diffusion  layer  and  the  catalyst  layer  can  be  deter¬ 
mined  in  terms  of  the  porosity  based  on  Bruggemann’s  equation. 
The  equations  are  listed  below: 


^CP  —  ^CP, 

(20) 

^GDL  =  (1  -  £GDL)0GDL  X  ^GDL. 

(21) 

acl  =  4ccl  x  o’s.cl. 

(22) 

Here,  £SiCl  and  <7SiCl  denote  the  volume  fraction  (porosity)  and 
electronic  conductivity,  respectively,  of  solid  catalyst  particles  in 
the  catalyst  layer.  In  addition,  the  source  term  §d  in  Eq.  (10)  for  the 
carbon  plate,  porous  gas  diffusion  layer  and  catalyst  layer  can  be 
expressed  as 


£cr  =  0,  (23) 

§GDL  =  0,  (24) 

(25> 

Boundary  conditions  at  the  anode  flow  and  the  cathode  flow 
channels  are  as  follows:  the  inlet  flow  rates  are  constant,  the  inlet 
gas  compositions  are  constant,  and  the  flows  are  fully  developed 
at  the  outlets  of  the  anode  and  cathode  flow  channels.  At  the  solid 
walls,  non-slip  and  zero  fluxes  hold.  At  the  interfaces  between  the 
gas  channels,  the  diffuser  layers,  the  catalyst  layers,  and  the  PEM, 
equality  is  applied  to  the  velocity,  mass  fraction,  momentum  flux, 
and  mass  flux.  In  this  study,  all  the  above  conservation  equations 
and  property  equations  are  solved  by  utilizing  a  general  purpose 
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computational  fluid  dynamic  code,  CFD-ACE+,  employing  a  finite 
volume  method,  released  by  ESI-CFD  Inc.  [19],  since  it  has  the  ability 
to  handle  this  moving  boundary  problem.  The  “moving  boundary 
problem”  actually  implies  that  the  boundary  of  the  gas  channel  is 
subject  to  change  during  each  iteration  of  the  process. 

The  coupled  set  of  equations  was  solved  iteratively,  and  the 
solution  was  considered  to  be  convergent  when  the  relative  error 
in  each  field  between  two  consecutive  iterations  was  less  than  a 
specified  small  number. 

3.  Shape  generation  for  gas  channel 


Inlet 


/ 


/ 


The  reason  for  using  a  B-spline  curve  for  optimal  channel  shape 
is  that  it  can  generalize  a  smooth  curve  that  satisfies  the  desired 
current  distribution  and  makes  it  a  proper  shape  for  a  channel; 
therefore  the  coordinates  of  the  control  points  become  the  param¬ 
eters  that  need  to  be  estimated  in  this  shape  design  problem  to 
determine  the  optimal  gas  channel  geometry. 

Let  us  consider  a  Cartesian  product  parametric  B-spline  curve 
[20]  given  by 


n+ 1 

nt) = hnin  —  £  —  fmaxi  2  <  k<n  + 1,  (26) 

j= i 


where  Bj  are  the  position  vectors  of  the  n  + 1  defining  polygon  ver¬ 
tices  and  Njfc(t)  are  the  normalized  B-spline  basis  functions.  For 
the  jth  normalized  B-spline  basis  function  of  order  k  (degree  k  —  1 ), 
the  basis  functions  Nj  k(t)  are  defined  by  the  Cox-deBoor  recursion 
formulas.  Specifically: 


and 


J  1  if  Xj  <t  <  Xj+ 1 

1  0  otherwise 

(*  -  (Xj+k  ~  t)Nj+hk-i(t) 

xj+k-\  ~  xj  xj+k  ~  xj+ 1 


(27) 


where  the  values  of  x2  are  the  elements  of  a  uniform  knot  vector 
satisfying  the  relation  Xj<t<  xJ+i .  The  parameter  t  varies  from  tmin 
to  tmax  along  the  curve  *P(t),  i.e.,  t )  are  the  curve  data  points. 

The  N  basis  functions  can  be  determined  from  the  knot  vector 
and  the  parameter  value  t.  Bj  are  the  required  polygon  net  points 
(control  points).  If  Bj  are  given,  the  curve  data  points  ^(t)  can  be 
calculated  from  Eq.  (26). 


4.  The  inverse  design  problem 

The  current  density  distributions  in  the  outlet  region  of  the 
gas  channel  are  almost  uniform  in  the  channel’s  width  direction, 
therefore  the  redesigned  gas  channel  surface  can  be  assumed  as 
a  two-dimensional  surface,  i.e.,  it  varies  only  in  the  length  (flow) 
and  height  directions,  and  the  dimensions  are  all  the  same  in  the 
channel’s  width  direction.  This  implies  that  we  need  to  redesign 
the  optimal  shape  of  the  gas  channel  only  in  the  length-height 
cross-section,  and  this  result  is  used  to  generate  the  redesigned 
gas  channel  surface  in  the  width  direction. 

For  the  inverse  geometry  design  problem,  the  shape  of  the 
redesigned  portion  of  the  gas  channel  is  regarded  as  being  unknown 
and  controlled  by  a  set  of  control  points;  in  addition,  the  desired 
distributions  of  current  density  located  on  the  upper  CP  surface  of 
the  redesigned  gas  channel  at  the  cathode  are  considered  available. 

Let  the  desired  current  densities  located  on  the  upper  CP  surface 
of  the  redesigned  gas  channel  at  the  cathode  be  denoted  by  Y(x2, 
yi)  =  Yj,  i  =  1  to  M,  where  M  represents  the  number  of  grid  points  for 
the  redesigned  portion  of  the  gas  channel.  The  inverse  geometry 
design  problem  can  then  be  stated  as  follows:  utilizing  the  above 


Outlet 


in 


Carbon  Plate  » 

Cathode  Gas  Channel  »  CGC 
Cathode  Gas  Diffusion  Layer  *  CGDL 
Cathode  Catalyst  Layer  »  CCL 
Proton  Exchange  Membrane  *  PEM 
Anode  Catalyst  Layer  »  ACL 
Anode  Gas  Diffusion  Layer  »  AGDL 
Anode  Gas  Channel  »  AGC 
Carbon  Plate  »  CP 


Redesigned  position  of 
the  gas  channel 


Fig.  1.  Physical  model  of  a  single-cell,  parallel  straight  channel  PEMFC. 


mentioned  desired  current  densities  Y2,  design  the  new  shape  for 
the  gas  channel. 

The  solution  of  the  present  inverse  geometry  design  problem 
is  to  be  obtained  in  such  a  way  that  the  following  functional  is 
minimized: 

M 

jmBj)}  =^[/c,i(Bj)  —  Yj]2  =  l/TU;  j  —  1  to  P.  (28) 

1=1 

Here,  /C2  represent  the  estimated  or  computed  current  densities 
at  the  channel  locations  (x2,  y2).  These  quantities  are  determined 
from  the  solution  of  the  direct  problem  given  previously  by  using 
an  estimated  gas  channel  ^(Bj),  with  P  representing  the  number  of 
control  points. 


4.1.  The  Levenberg-Marquardt  method  (LMM)  for  minimization 


If  the  redesigned  portion  of  the  gas  channel  is  discretized  into 
M  points  and  P  control  points  are  used,  Eq.  (28)  is  minimized  with 
respect  to  the  estimated  parameters  Bj  to  obtain: 


mnBj)] 

dBj 


E 


3/c  ,t 
dBj 


Yi]  =0:  j  =  l  toP, 


(29) 


where  M  should  be  equal  to  or  greater  than  P,  otherwise  an 
underdetermined  system  of  equations  will  be  obtained,  and  it 
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Table  1 

Fixed  parameters  used  by  CFD-ACE+  code  in  the  present  study. 


Parameters 

Symbols 

Values 

Reference  transfer  current  density  at  anode  **(Am-2) 

4, ref 

9E+008 

Reference  transfer  current  density  at  cathode  (A  m-2 ) 

4, ref 

2.5E+002 

Anodic  charge  transfer  coefficients  for  anode  reaction 

Ota, a 

0.5 

Cathodic  charge  transfer  coefficients  for  anode  reaction 

Ota,c 

0.5 

Anodic  charge  transfer  coefficients  for  cathode  reaction 

Otc,a 

1.5 

Cathodic  charge  transfer  coefficients  for  cathode  reaction 

Otc,c 

1.5 

Thickness  of  carbon  plate  (m) 

dCp 

2E-003 

Width  of  rib  (m) 

Wr 

IE-003 

Resistivity  of  carbon  plate  (ohm-m) 

J2cp 

2.7E-004 

Length  of  gas  channel  (m) 

Leu 

5E-002 

Height  of  gas  channel  (m) 

Lieu 

IE-003 

Width  of  gas  channel  (m) 

WCu 

IE— 003 

Anode  inlet  gas  velocity  (m  s-1 ) 

Va 

0.3 

Cathode  inlet  gas  velocity  (ms-1) 

VC 

0.5 

Anode  inlet  mass  fraction  of  H2 

ZU2,a 

0.445 

Anode  inlet  mass  fraction  of  H20 

ZU20,a 

0.555 

Cathode  inlet  mass  fraction  of  02 

Z0  2,c 

0.212 

Cathode  inlet  mass  fraction  of  N2 

ZN  2,c 

0.709 

Cathode  inlet  mass  fraction  of  H20 

ZU20,c 

0.079 

Thickness  of  gas  diffusion  layer  (m) 

dGDL 

0.3E-003 

Porosity  of  gas  diffusion  layer  (m2) 

£gdl 

0.5 

Permeability  of  gas  diffusion  layer  (m2) 

/<GDL 

1.76E— 011 

Bruggemann  coefficient  of  gas  diffusion  layer 

#GDL 

1.5 

Electronic  conductivity  of  gas  diffusion  layer  (£2-1  m-1 ) 

^GDL 

300 

Thickness  of  catalyst  layer  (m) 

dei 

IE— 005 

Porosity  of  catalyst  layer  (m2) 

£ei 

0.112 

Permeability  of  catalyst  layer  (m2) 

kc  L 

1.76E— 011 

Bruggemann  coefficient  of  catalyst  layer 

0CL 

1.5 

Electronic  conductivity  of  catalyst  layer  (£2-1  m-1 ) 

tfCL 

300 

Thickness  of  membrane  (m) 

dm 

5E-005 

Porosity  of  membrane  (m2) 

£m 

0.28 

Permeability  of  membrane  (m2) 

km 

5E-018 

Bruggemann  coefficient  of  membrane 

13 

will  be  impossible  to  calculate  the  inverse  solutions.  Eq.  (29) 
is  linearized  by  expanding  Ic,i{Bj)  in  Taylor  series  and  retaining 
the  first-order  terms.  Then  a  damping  parameter  /xn  is  added  to 
the  resulting  expression  to  improve  convergence,  leading  to  the 
Levenberg-Marquardt  method  [21  ],  given  by 

(F  +  /znE)AB  =  D  (30) 


Fig.  3.  The  original  and  estimated  optimal  shapes  of  gas  channels  at  (a)  0.4  V  and 
(b)  0.7  V,  with  20%  increase  in  the  local  current  density. 


F  =  ftTft 


(31)  D  =  ftTU 


(32) 


AB  =  Bn+1-B".  (33) 

Here,  the  superscripts  n  and  T  represent  the  iteration  index  and 
transpose  matrix,  respectively,  E  is  the  identity  matrix,  and  ft 
denotes  the  Jacobian  matrix,  defined  as 

(34) 

9Bt 

The  Jacobian  matrix  defined  by  Eq.  (34)  is  determined  by  per¬ 
turbing  the  unknown  parameters  Bj  one  at  a  time  and  computing 
the  resulting  change  in  current  density  from  the  solution  of  the 
direct  problem,  Eqs.  (1  )-(25). 

Eq.  (30)  is  now  written  in  a  form  suitable  for  iterative  calculation 
as 

Bn+1  =B"  +  (dT'»  +  AtnSr1'dT(Ic-Y).  (35) 

When  /in  =  0,  Newton’s  method  is  obtained,  and  as  /xn  -*  oo, 
the  steepest-descent  method  is  obtained.  For  fast  convergence  the 
steepest-descent  method  is  applied  first,  then  the  value  of  /xn 
is  decreased,  and  finally  Newton’s  method  is  used  to  obtain  the 
inverse  solution.  The  algorithm  for  choosing  this  damping  value  /xn 
is  described  in  detail  by  Marquardt  [21  ],  so  it  is  not  repeated  here. 

The  bridge  between  CFD-ACE+  and  LMM  is  the  INPUT/OUTPUT 
files.  These  files  should  be  arranged  such  that  their  format  can  be 
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0.05  -n 

H  (b)  0.05  - 1 

(c)  0.05  — 

"1  (d)  0.05 -T 

0.045- 

0.045- 

0.045 -i 

0.045- 

0.04- 

0.04  - 

0.04 -j 

0.04- 

0.035- 

0.035- 

0.035  — 

0.035- 

0.03- 

m 

0.03- J 

0.03- 

0.03- 

0.025- 

0.025- 

0.025- 

0.025- 

0.02- 

0.02- 

0.02- 

0.02- 

0.015-  ( 

0.015- 

0.015-1 

0.015- 

0.01- 

0.01- 

i— i 

o 

o 

o.oi- 

0.005  - 

0.005- 

0.005- 

0.005- 

Y  0-1 

ux 

®  o-l 

1  0_l 

i  0- 

Jsoliclz  -  A/m2 
1.484E+004 


6000- 


5975 


Fig.  4.  The  current  density  distributions  for  (a)  original  case,  (b)  Case  A,  (c)  Case  B  and  (d)  Case  C  gas  channels  at  0.4  V,  with  20%  increase  in  the  local  current  density. 


(a)  0.05- 

0.045- 

0.04- 

0.035- 

0.03- 
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0.02- 

0.015- 

0.01- 
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0.025- 
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0.015- 
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(d) 

0.035- 

0.03- 

0.025- 

0.02- 

0.015- 
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Jsolicl_z-A/m2 

4350 


3266 


Fig.  5.  The  current  density  distributions  for  (a)  original  case,  (b)  Case  A,  (c)  Case  B  and  (d)  Case  C  gas  channels  at  0.7  V,  with  20%  increase  in  the  local  current  density. 
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Length  (m) 

Fig.  6.  The  original  and  calculated  (a)  saturations  and  (b)  oxygen  concentrations 
along  the  center  line  of  the  cathode  GDL-CL  interface  at  0.4  V,  with  20%  increase  in 
the  local  current  density. 


recognized  by  CFD-ACE+  and  LMM.  A  sequence  of  forward  PEMFC 
problems  is  solved  by  CFD-ACE+  in  an  effort  to  update  the  gas  chan¬ 
nel  geometry  by  minimizing  a  residual  measuring  the  difference 
between  estimated  and  desired  current  densities  located  on  the 
upper  CP  surface  of  the  redesigned  gas  channel  at  the  cathode  under 
the  present  algorithm. 


5.  Computational  procedure 

The  iterative  computational  procedure  for  the  solution  of  this 
inverse  design  problem  using  the  Levenberg-Marquardt  method 
can  be  summarized  as  follows: 

Step  1.  Choose  the  initial  guess  for  control  points  B  (obtained  by 
using  the  original  shape  of  the  gas  channel  and  B-spline  curve 
fitting)  at  iteration  n  to  start  the  computations. 

Step  2.  Solve  the  direct  problem  given  by  Eqs.  (l)-(25)  to  obtain 
computed  current  densities  lc. 

Step  3.  Construct  the  Jacobian  matrix  in  accordance  with  Eq.  (34). 
Step  4.  Update  B  from  Eq.  (35). 

Step  5.  Check  the  stopping  criterion;  if  not  satisfied  go  to  Step  2 
and  iterate. 


Length  (m) 


Length  (m) 

Fig.  7.  The  original  and  calculated  (a)  saturations  and  (b)  oxygen  concentrations 
along  the  center  line  of  the  cathode  GDL-CL  interface  at  0.7  V,  with  20%  increase  in 
the  local  current  density. 


6.  Results  and  discussion 

The  accuracy  of  the  solution  from  CFD-ACE+  plays  a  very 
important  role  in  this  inverse  design  problem.  If  the  solution 
from  CFD-ACE+  cannot  reproduce  the  actual  performance  of  the 
PEMFC,  the  inverse  solutions  can  never  be  obtained  accurately. 
The  first  task  is  thus  to  show  the  validity  of  the  solutions 
from  CFD-ACE+  by  comparing  them  with  existing  experimental 
data. 

The  benchmark  problem  for  the  numerical  solution  using  CFD- 
ACE+  is  examined  here  based  on  the  experiments  conducted  by 
Cheng  et  al.  [16].  The  grid  system  for  a  single-cell  fuel  cell  module 
with  a  cross-sectional  area  of  22  cm  x  22  cm  and  an  active  area  of 
14  cm  x  14  cm  has  been  constructed.  Two  2-mm  thick  composite 
graphite  plates  are  used  as  the  flow  fields  and  current  conductors. 
A  five-layer  membrane  exchange  assembly  (MEA)  is  placed  at  the 
center  between  the  two  composite  graphite  plates.  Fig.  1  shows 
the  present  single-cell  PEMFC.  For  simplicity,  only  one  gas  chan¬ 
nel  is  computed  in  the  present  study,  as  in  Cheng  et  al.  [16].  The 
inlet  velocities  for  the  hydrogen  and  the  air  are  chosen  as  0.3  and 
0.5  m  s_1,  respectively.  The  hydrogen  and  the  air  are  humidified  to 
reach  100%  relative  humidity  and  80  °C  temperature  before  enter¬ 
ing  the  gas  flow  channels.  The  operating  and  geometrical  conditions 
given  to  the  direct  problem  solver  are  consistent  with  the  experi¬ 
ments.  The  fuel  cell  temperature  is  maintained  at  70  °C.  The  fixed 
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Fig.  8.  The  pressure  distributions  for  (a)  original  case,  (b)  Case  A,  (c)  Case  B  and  (d)  Case  C  gas  channels  at  0.4  V,  with  20%  increase  in  the  local  current  density. 


parameters  used  for  the  numerical  calculations  in  CFD-ACE+  are 
listed  in  Table  1. 

Fig.  2  illustrates  the  comparisons  between  the  numerical  solu¬ 
tions  using  CFD-ACE+  and  the  experimental  data  from  [16]  for 
the  polarization  curves.  It  is  obvious  that  the  numerical  solutions 
from  CFD-ACE+  match  well  with  the  experimental  data.  Three  gas 
channel  design  problems  will  be  considered  here,  differing  by  the 
starting  point  of  the  redesigned  channel: 

(1)  Case  A:  the  redesigned  channel  length  equals  20  mm  (from  30 
to  50  mm), 

(2)  Case  B:  the  redesigned  channel  length  equals  15  mm  (from  35 
to  50  mm),  and 


(3)  Case  C:  the  redesigned  channel  length  equals  10  mm  (from  40 
to  50  mm). 

The  number  of  control  points  is  taken  as  four  in  all  the  examples 
considered  here,  and  the  design  criterion  is  that  the  current  density 
on  the  upper  CP  surface  of  the  redesigned  channel  is  required  to 
increase  by  20%  of  the  original  current  density  (for  a  rectangular  gas 
channel)  at  0.4  V  (low  operating  voltage)  and  0.7  V  (high  operating 
voltage).  One  should  note  that  even  though  more  control  points 
describe  the  unknown  surface  more  accurately,  they  also  require 
more  computer  time  to  obtain  the  inverse  solutions. 

The  inverse  calculations  are  performed  by  using  the 
Levenberg-Marquardt  method  based  on  the  above  conditions.  The 
initial  estimates  of  Bj  are  obtained  via  B-spline  curve  fitting  for 


Fig.  9.  The  pressure  distributions  for  (a)  original  case,  (b)  Case  A,  (c)  Case  B  and  (d)  Case  C  gas  channels  at  0.7  V,  with  20%  increase  in  the  local  current  density. 


Voltage  (V)  Original  total  current  Case  A  (0.4  V)  total  current  Case  B  (0.4  V)  total  current  Case  C  (0.4  V)  total  current  Case  A  (0.7  V)  total  current  Case  B  (0.7  V)  total  current  Case  C  (0.7  V)  total  current 

density  (Am-2)  density  (Am-2)  density  (Am-2)  density  (Am-2)  density  (Am-2)  density  (Am-2)  density  (Am-2) 
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(a) 


Fig.  10.  The  original  and  estimated  (a)  optimal  shape  of  gas  channels  and  (b)  current 
density  distributions  at  0.4  V,  with  30%  increase  in  the  local  current  density. 


the  rectangular  gas  channel.  After  only  7-10  iterations  (the  CPU 
time  on  an  Intel  Dual-core  1.8  GHz  processor  is  about  3-4  days), 
the  estimated  shape  of  the  gas  channel  can  be  obtained,  and  this 
result  is  shown  in  Fig.  3(a)  and  (b).  The  average  relative  errors  for 
the  desired  and  calculated  current  density  on  the  upper  surface 
of  the  redesigned  channel  are  calculated  as  2.5%  (Case  A  at  0.4  V), 
2.4%  (Case  B  at  0.4  V),  2.5%  (Case  C  at  0.4  V),  2.7%  (Case  A  at  0.7  V), 
2.7%  (Case  B  at  0.7  V),  and  2.6%  (Case  C  at  0.7  V),  where  the  average 
relative  error  ERR  is  defined  as 


M 


ERR  =  £ 


hi  ~  Yi 


100 -M, 


(36) 


where  M  =  120, 90  and  60  for  cases  A,  B  and  C,  respectively.  It  is  obvi¬ 
ous  from  these  errors  that  the  Levenberg-Marquardt  method  has 
been  applied  successfully  for  estimating  the  optimal  shape  of  gas 
channels  in  the  numerical  examples,  since  the  calculated  current 
densities  closely  agree  with  the  desired  ones.  The  overall  increased 
percentages  for  the  current  in  these  six  cases  are  obtained  as  3.9% 
(Case  A  at  0.4  V),  2.5%  (Case  B  at  0.4  V),  1.6%  (Case  C  at  0.4  V),  4.8% 
(Case  A  at  0.7  V),  3.6%  (Case  B  at  0.7  V),  and  2.4%  (Case  C  at  0.7  V), 
respectively. 

Fig.  3(a)  and  (b)  shows  the  original  and  redesigned  optimal 
shapes  for  gas  channels  at  operating  voltages  of  0.4  and  0.7  V, 
respectively.  The  redesigned  optimal  gas  channels  exhibit  oscilla¬ 
tory  profiles,  and  the  first  point  of  the  redesigned  channel  is  always 
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Fig.  11.  The  original  and  calculated  (a)  saturations  and  (b)  oxygen  concentrations 
along  the  center  line  of  the  cathode  GDL-CL  interface  at  0.4  V,  with  30%  increase  in 
the  local  current  density. 


at  the  lowest  position.  This  may  be  due  to  the  fact  that  at  the  first 
point,  a  sudden  increase  in  the  desired  current  density  is  required; 
to  obey  this  requirement,  a  sudden  contraction  of  the  gas  channel  is 
needed.  However,  at  the  adjacent  section  of  the  channel,  the  height 
does  not  have  to  be  corrected  as  low  as  that  in  the  first  position  to 
match  the  desired  current  density;  therefore,  the  profile  of  the  gas 
channel  increases  in  height  smoothly,  and  then  decreases  in  height 
smoothly  if  necessary. 

To  show  that  the  estimated  shape  is  indeed  the  optimal  shape, 
we  have  considered  an  artificial  rectangular  gas  channel  with  a 
sudden  contraction  and  the  same  duct  volume  as  Case  C  at  0.4  V. 
Fig.  3(a)  shows  the  profile  for  this  artificially  contracted  gas  chan¬ 
nel.  The  direct  solution  is  then  calculated  based  on  this  gas  channel, 
and  the  average  relative  error  for  the  desired  and  calculated  current 
density  on  the  upper  surface  of  the  redesigned  channel  is  calculated 
as  2.7%,  which  is  greater  than  that  for  the  optimal  gas  channel.  The 
overall  increase  in  the  current  is  1.4%,  which  is  also  smaller  than  that 
for  the  optimal  gas  channel.  It  can  therefore  be  concluded  that  the 
performance  of  the  PEMFC  with  an  optimal  redesigned  gas  channel 
is  indeed  better  than  its  performance  with  the  artificial  gas  channel. 
Figs.  4  and  5  indicate  the  original  and  calculated  current  density 
on  the  cathode  CP  surface  for  Cases  A,  B  and  C  at  0.4  and  0.7  V, 
respectively.  It  is  obvious  from  these  two  figures  that  the  current 
densities  in  the  redesigned  gas  channels  are  all  higher  than  those 
for  the  original  gas  channel. 

From  Fig.  3(a)  and  (b),  we  found  that  in  order  to  match  the 
desired  current  density,  the  gas  channel  must  be  contracted  sud¬ 
denly  at  the  start  of  the  redesigned  position  of  the  gas  channel, 
in  order  to  remove  the  saturated  water  as  well  as  to  increase  the 
concentration  of  oxygen.  The  above  statements  can  be  verified  by 
plotting  the  calculated  water  saturation  and  oxygen  concentration 
with  channel  length.  Here,  water  saturation  is  defined  as  the  ration 
of  the  liquid  water  volume  in  the  porous  medium  to  the  total  vol¬ 
ume  of  the  porous  medium.  Fig.  6(a)  and  (b)  indicate  the  original 
and  calculated  saturation  and  oxygen  concentration  along  the  cen¬ 
ter  line  of  the  gas  channel  and  on  the  cathode  GDL-CL  interface, 
respectively,  for  Cases  A,  B  and  C,  at  the  low  operating  voltage  of 
0.4  V,  while  Fig.  7(a)  and  (b)  shows  the  same  for  the  high  operating 
voltage  of  0.7  V. 


UY  Length  (m) 


Fig.  12.  The  pressure  distributions  for  (a)  original  case,  (b)  Case  A,  (c)  Case  B  and  (d)  Case  C  gas  channels  at  0.4  V,  with  30%  increase  in  the  local  current  density. 
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Table  3 

Polarization  data  for  the  original  and  various  redesigned  channels  with  30%  increase  in  the  local  current  density. 


Voltage,  V 

Original  total  current 
density  (Am-2) 

Case  A  (0.4  V)  total 
current  density  (A  m-2 ) 

Case  B  (0.4  V)  total 
current  density  (A  m-2 ) 

Case  C  (0.4  V)  total 
current  density  (A  m-2 ) 

0.97 

285 

456 

416 

386 

0.95 

404 

526 

484 

452 

0.92 

562 

696 

650 

615 

0.89 

785 

932 

882 

843 

0.85 

1158 

1,321 

1,265 

1,222 

0.83 

1383 

1,555 

1,496 

1,451 

0.8 

1785 

1,967 

1,905 

1,857 

0.75 

2653 

2,853 

2,785 

2,732 

0.7 

3775 

3,996 

3,921 

3,864 

0.65 

5094 

5,341 

5,258 

5,196 

0.6 

6494 

6,785 

6,687 

6,617 

0.575 

7171 

7,499 

7,393 

7,313 

0.55 

7781 

8,164 

8,045 

7,956 

0.5 

8677 

9,190 

9,032 

8,909 

0.4 

9482 

10,115 

9,917 

9,760 

0.3 

9729 

10,381 

10,156 

10,005 

From  Figs.  6(a)  and  7(a),  we  found  that  the  longer  the  redesigned 
channel,  the  more  saturated  water  will  be  removed,  and  that  the 
saturation  at  0.4  V  is  more  than  that  at  0.7  V.  Contrarily,  the  oxygen 
concentration  at  0.4  V  is  less  than  that  at  0.7  V.  This  is  because  the 
electrochemical  reaction  is  more  moderate  at  0.7  V  than  at  0.4  V, 
and  as  a  result  less  oxygen  will  be  consumed  and  less  water  will  be 
produced  at  0.7  V.  This  also  explains  why  the  overall  increase  of  cur¬ 
rent  at  0.4  V  is  higher  than  that  at  0.7  V.  From  Figs.  6(b)  and  7(b),  at 
the  redesigned  position  of  the  original  gas  channel,  the  oxygen  con¬ 
centration  for  0.4  V  is  much  lower  than  that  for  0.7  V.  This  implies 
a  perfect  reaction  at  0.4  V.  By  imposing  the  requirement  that  the 
redesigned  current  density  be  20%  greater  than  the  original  current 
density  at  the  redesigned  position,  the  amount  of  the  increment  of 
oxygen  concentration  at  0.4  V  is  also  almost  consumed  when  com¬ 
pared  with  the  case  at  0.7  V,  and  therefore  the  overall  increase  of 
current  at  0.4  V  is  greater  than  that  at  0.7  V. 

The  pressure  is  also  an  important  index  for  this  design  problem, 
since  it  will  be  increased  for  the  optimal  gas  channel.  Figs.  8  and  9 
indicate  the  pressure  distributions  along  the  gas  channel  at  0.4  and 
0.7  V,  respectively,  for  the  original  case,  Case  A,  Case  B  and  Case  C  at 
0.4  and  0.7  V.  It  is  obvious  that  the  longer  the  redesigned  gas  chan¬ 
nel,  the  larger  the  pressure  difference.  If  this  pressure  drop  cannot 
be  handled  by  the  gas  pumping  system,  the  optimally  designed  gas 
channel  cannot  be  used. 

Table  2  shows  the  polarization  data  based  on  the  optimal  gas 
channels,  Case  A,  Case  B  and  Case  C,  designed  for  0.4  and  0.7  V,  when 
operating  at  any  other  voltages.  It  is  obvious  that  these  designs 
are  for  0.4  and  0.7  V,  but  can  still  be  applied  to  different  operating 
voltages  to  ensure  increase  of  the  overall  current. 

The  next  question  is  whether  the  current  density  on  the  upper 
CP  surface  of  the  redesigned  channel  can  be  required  to  increase  by 
30%  of  the  original  current  density  (for  a  rectangular  gas  channel) 
at  0.4  V  (low  operating  voltage)  and  0.7  V  (high  operating  voltage). 
The  inverse  calculations  are  performed  again,  using  the  Levenberg- 
Marquardt  method  based  on  the  above  conditions. 

Unfortunately,  at  0.7  V,  the  design  fails,  since  it  is  impossible 
to  obtain  the  needed  shape  to  match  the  desired  current  density. 
However,  at  0.4  V,  it  can  be  estimated  successfully.  Fig.  10(a)  shows 
the  original  and  estimated  optimal  shapes  for  gas  channel  at  0.4  V. 
The  tendency  of  the  designed  optimal  gas  channels  is  similar  to  the 
previous  estimations.  The  average  relative  errors  for  the  desired 
and  calculated  current  densities  on  the  upper  CP  surface  of  the 
redesigned  channel  are  calculated  as  2.8%  (Case  A  at  0.4  V),  2.6% 
(Case  B  at  0.4  V)  and  2.7%  (Case  C  at  0.4  V),  and  again  the  calculated 
current  densities  closely  agree  with  the  desired  densities.  The  over¬ 
all  increased  percentages  for  the  current  in  the  present  cases  are 


obtained  as  6.7%  (Case  A  at  0.4  V),  4.6%  (Case  B  at  0.4  V)  and  2.9% 
(Case  Cat  0.4  V). 

Fig.  10(b)  indicates  the  original  and  calculated  current  densities 
for  Cases  A,  B  and  C  at  0.4  V.  Fig.  11(a)  and  (b)  shows  the  original  and 
calculated  saturation  and  oxygen  concentration  along  the  center 
line  of  the  gas  channel  and  on  the  cathode  GDL-CL  interface  for 
Cases  A,  B  and  C  at  0.4  V.  By  comparing  Fig.  11(a)  with  Fig.  7(a)  and 
Fig.  11(b)  with  Fig.  7(b),  it  can  be  seen  that  more  water  is  removed 
and  higher  oxygen  concentration  is  obtained  under  this  redesigned 
condition.  Fig.  12  indicates  the  pressure  distributions  along  the  gas 
channel  for  the  original  case,  Case  A,  Case  B  and  Case  C  at  0.4  V.  It 
is  obvious  by  comparing  Figs.  8  and  12  that  the  higher  the  desired 
current  density,  the  larger  the  pressure  difference. 

Table  3  indicates  the  polarization  data  based  on  the  optimal 
gas  channels,  Case  A,  Case  B  and  Case  C,  designed  at  0.4  V,  when 
operating  at  any  other  voltages.  Again,  the  design  at  0.4  V  can 
still  be  applied  to  different  operating  voltages  to  obtain  better 
performance  from  the  fuel  cell.  Based  on  the  above  numerical 
experiment,  it  is  concluded  that  the  design  problem  of  determin¬ 
ing  the  optimal  shape  for  a  PEMFC  fuel  cell  gas  channel  using  the 
Levenberg-Marquardt  method  is  now  completed.  The  calculated 
current  density  can  match  well  with  the  desired  value  if  the  require¬ 
ment  is  attainable. 

7.  Conclusions 

An  inverse  design  problem  in  estimating  the  optimal  shape 
of  gas  channels  for  PEMFCs  from  the  knowledge  of  desired  cur¬ 
rent  density  distributions,  using  the  techniques  of  B-spline  curves 
and  the  Levenberg-Marquardt  method,  has  been  developed  and 
applied  successfully.  Numerical  experiments  have  been  performed 
based  on  low  and  high  operating  voltages,  in  this  case  0.4  and  0.7  V. 
The  results  reveal  that  the  present  algorithm  needs  only  a  few  itera¬ 
tions  to  obtain  the  optimal  shape  of  the  gas  channel,  and  the  overall 
current  can  be  increased  by  using  the  estimated  gas  channel. 

The  advantages  of  using  the  technique  of  the  inverse  design 
problem  in  designing  the  optimal  shape  of  gas  channels  lie  in  the 
facts  that  (1 )  the  time  needed  to  redesign  the  channel  can  be  short¬ 
ened  and  (2)  the  desired  current  density  can  be  well  matched  when 
compared  with  the  traditional  trial  and  error  method. 
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